rm(list = ls(all.names = TRUE))
gc()
set.seed(211917)
options(scipen=999)

packages <-c("tidyverse","foreign","stargazer","naniar","ri2")

new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

lapply(packages, require, character.only = TRUE)
rm(packages, new.packages)

setwd("PUT YOUR WORKING DIRECTORY HERE")

# load data
dat <- read.dta("./data.dta")

# some preliminary cleaning
dat$q6_years_married[dat$q6_years_married == -9] <- 0
dat$treatment_binary <- ifelse(dat$holy_endorser==1 | dat$secular_endorser==1, 1, 0)
dat$q30_combined_s <- scale(dat$q30_combined)

#---------------#
# Data cleaning #
#---------------#

## select all possible covariates and outcome at first
dat <- dat %>% 
  select(q30_combined, q30_combined_s, treatment_binary, respondent_muslim, respondent_male, enumerator_muslim, enumerator_male, q1_age, q2_male, q3_kamba, q3_kikuyu, q3_luo, q3_somali, q3_other, q4_religion_christian, q4_religion_muslim, q5_convert, q6_years_married, q6_married, q7_has_job, q8_income_last_month, q8_has_income, q9_chance_become_rich, q10_plans_to_vote, q11_govt_represents_interest, q12_feels_on_loosing_side, q13_a_knows_emigrant, q13_b_knows_emigrant_som_eri, q15_witnesses_interrel_violence, q16_witnesses_viol_musl_govt, q17_friend_family_died, q18_lost_job, q18_arrested, q18_house_of_worship_raided, q18_stopped_speaking_parents, q18_love_problems, q18_friend_left_country, q18_sum, q19_rlts_mother, q19_not_raised_by_mom, q19_mother_dead, q20_rlts_father, q20_not_raised_by_father, q20_father_dead, q21_respect_friends_family, q22_gender, q22_religion, q22_tribal, q22_national, q22_youth, q23_religious_attendance, q23_frequents_house_worship, q24_deontology, q25_payment_acceptable, q26_payment_ethical, q27_radicalization, q33_distance)

## check for missings
vis_miss(dat)

## following variables seems problematic, i.e., shouldn't impute at this level
table(is.na(dat$q13_a_knows_emigrant))
table(is.na(dat$q19_rlts_mother))
table(is.na(dat$q20_rlts_father))
table(is.na(dat$q23_religious_attendance))

## kick problematic vars out
dat <- dat %>% select(-q13_a_knows_emigrant, -q19_rlts_mother, -q20_rlts_father)

## check missing again:
vis_miss(dat) 

## remove one missing treatment indicator (can't be imputed)
dat <- dat[which(dat$treatment_binary>-1),]

## impute rest of missings with mean 
for(i in 1:ncol(dat)){
  dat[is.na(dat[,i]), i] <- mean(dat[,i], na.rm = TRUE)
}
rm(i)

## 
dat <- dat[complete.cases(dat), ]

#-------------------------------------------------#
# OLS with all covariates and FEs and interations #
#-------------------------------------------------#

# set seed
set.seed(12345)

## extract all possible controls
controls <- as.character(colnames(dat[,3:54]))

## religion
dat$it_mus <- dat$treatment_binary * dat$respondent_muslim
lm_muslim <- lm(paste("q30_combined ~ it_mus + ", paste(controls, collapse=" + ")), data = dat)
lm_muslim

# any attendance
dat$it_frhw <- dat$treatment_binary * dat$q23_frequents_house_worship
lm_attendance <- lm(paste("q30_combined ~ it_frhw + ", paste(controls, collapse=" + ")), data = dat)

# identification
dat$it_relig <- dat$treatment_binary * dat$q22_religion
lm_identity <- lm(paste("q30_combined ~ it_relig + ", paste(controls, collapse=" + ")), data = dat)

# convert
dat$it_conv <- dat$treatment_binary * dat$q5_convert
lm_convert <- lm(paste("q30_combined ~ it_conv + ", paste(controls, collapse=" + ")), data = dat)

# social distance
dat$it_dist <- dat$treatment_binary * dat$q33_distance
lm_distance <- lm(paste("q30_combined ~ it_dist + ", paste(controls, collapse=" + ")), data = dat)

# extremism
dat$it_rad <- dat$treatment_binary * dat$q27_radicalization
lm_extremism <- lm(paste("q30_combined ~ it_rad + ", paste(controls, collapse=" + ")), data = dat)

# witnessed violence
dat$it_intrelvio <- dat$treatment_binary * dat$q15_witnesses_interrel_violence
lm_violence_inter <- lm(paste("q30_combined ~ it_intrelvio + ", paste(controls, collapse=" + ")), data = dat)

dat$it_viomusgov <- dat$treatment_binary * dat$q16_witnesses_viol_musl_govt
lm_violence_musl <- lm(paste("q30_combined ~ it_viomusgov + ", paste(controls, collapse=" + ")), data = dat)

dat$it_raid <- dat$treatment_binary * dat$q18_house_of_worship_raided
lm_violence_raid <- lm(paste("q30_combined ~ it_raid + ", paste(controls, collapse=" + ")), data = dat)


stargazer(lm_muslim, lm_attendance, lm_identity, lm_convert, lm_distance, lm_extremism, lm_violence_inter, lm_violence_musl, lm_violence_raid,
          keep = c("it_mus","it_frhw","it_relig","it_conv","it_dist","it_rad","it_intrelvio","it_viomusgov","it_raid"),
          covariate.labels = c("T * Muslim","T * Attends Church/Mosque","T * Religious identification","T * Convert","T * Social distance",
                               "T * Extremism","T * Interreligious violence","T * Govt. violence", "T * Church/Mosque raid"),
          add.lines = list(c("Controls", "Yes", "Yes", "Yes","Yes", "Yes", "Yes","Yes", "Yes", "Yes"),
                           c("FEs", "Yes", "Yes", "Yes","Yes", "Yes", "Yes","Yes", "Yes", "Yes")),
          keep.stat = c("n", "adj.rsq"),
          digits = 3,
          dep.var.caption = "Support for violence",
          dep.var.labels = c("","","","","","","","",""),
          label = "tab:main",
          title = "Effect heterogeneity",
          out = "PUT YOUR FILE PATH HERE")

